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Abstract — The  Naval  Oceanographic  Office  uses  the  Navy 
Coupled  Ocean  Data  Assimilation  (NCODA)  system  to  per¬ 
form  data  assimilation  for  ocean  modeling.  Currently  the 
system  uses  a  3D  multivariate  optimum  interpolation  (3D 
MVOI)  algorithm  to  produce  outputs  of  temperature,  salinity, 
geopotential,  and  u/v  velocity.  NCODA  is  run  in  a  stand¬ 
alone  mode  to  support  automated  ocean  data  quality  control 
(NCODA  OcnQC)  and  to  test  software  updates.  NCODA  is 
also  coupled  with  the  Regional/Global  Navy  Coastal  Ocean 
Model  (RNCOM/GNCOM).  The  RNCOM/NCODA  system 
is  being  used  as  part  of  an  Adaptive  Sampling  and  Predic¬ 
tion  (ASAP)  pre-operational  project,  that  makes  use  of  the 
Ensemble  Transform  (ET)  and  Ensemble  Transform  Kalman 
Filter  (ET  KF)  applied  to  ensemble  runs  of  the  RNCOM. 
The  ET  KF  is  used  to  predict  the  posterior  error  covariances 
resulting  from  possible  profile  measurements.  These  results 
aid  in  predicting  the  impact  of  ocean  observations  on  the 
future  analysis,  and  thus  allow  the  direction  of  limited  assets 
to  areas  that  will  have  the  maximum  gain  (for  applications 
such  as  ocean  acoustics).  A  review  of  these  systems  will 
be  given  as  well  as  examples  of  the  metrics  used  for  the 
RNCOM/NCODA  system,  ensemble  modeling,  and  ASAP. 

I.  NCODA 

The  Naval  Oceanographic  Office  uses  the  Navy  Cou¬ 
pled  Ocean  Data  Assimilation  (NCODA)  system  [1]  to 
perform  data  assimilation  for  ocean  modeling.  This  sys¬ 
tem  is  implemented  on  the  supercomputer  configura¬ 
tion  of  the  Navy  DoD  Supercomputing  Resource  Center 
(Navy  DSRC).  The  NCODA  analysis  is  running  in  both  a 
(1)  stand-alone  mode  (NCODA-SA)  to  support  real-time 
automated  ocean  data  quality  control  (NCODA  OcnQC) 
and  to  test  software  updates,  as  well  as  (2)  coupled 
with  the  Global  Navy  Coastal  Ocean  Model  (GNCOM) 
(GNCOM/NCODA)  and  Regional  Navy  Coastal  Ocean 
Model  (RNCOM)  (RNCOM/NCODA). 

The  NCODA  system  consists  of  a  frontend  ocean  data 
quality  control  (NCODA  OcnQC)  program  that  feeds 
an  analysis  program  (NCODA  Analysis).  Currently  the 


analysis  system  is  based  upon  a  3D  multivariate  opti¬ 
mum  interpolation  (3D  MVOI)  algorithm  [2]  to  produce 
3D  output  fields  of  temperature,  salinity,  geopotential, 
and  geostrophic  u/v  velocity  (T,  S,  <J>,  U,  V).  NCODA 
can  be  used  for  global  or  regional  applications  and  is 
relocatable  with  the  ability  to  have  multi-scale  analysis 
nests  (with  successively  higher  resolution  grids  in  a 
3:1  nest  ratio). 

A.  NCODA  Implementation 

The  NCODA  OcnQC  is  run  on  a  6-hour  schedule 
for  feeding  the  NCODA-SA,  RNCOM/NCODA,  and 
GNCOM/NCODA  systems.  Part  of  NCODA  OcnQC  is 
the  running  of  2D  OI  analyses  of  sea-surface  tempera¬ 
ture  (SST)  and  sea  ice  in  order  to  QC  SST  and  sea  ice 
observations.  These  are  performed  in  the  arctic  on  a 
9  km  polar  stereographic  grid  and  globally  on  a  9-12 
km  Mercator  grid  (midlatitude  versus  equatorial  grid 
spacing).  The  output  of  a  previous  NCODA  Analysis 
can  also  be  used  as  a  background  field  for  QC  of  in  situ 
profile  observations. 

The  NCODA  Analysis  (with  and  without  coupling 
with  an  ocean  model)  is  run  on  a  24  hour  schedule.  The 
global  analysis  is  performed  on  a  12-18  km  Mercator 
grid  (1/6°),  with  future  plans  of  9-13.5  km  grid  (1/8°). 
Regional  NCODAs  can  be  run  on  a  variety  of  grids,  but 
in  general  are  either  on  a  Mercator  or  spherical  projection 
at  a  nominal  3  km  grid  resolution. 

The  data  and  process  flow  of  the  NCODA  system  can 
be  seen  in  figure  1.  The  raw  observations  are  processed 
by  the  NCODA  OcnQC  and  assigned  "probability  of 
error"  (POE)  values.  The  data  and  its  associated  POE 
values  input  to  the  NCODA  Analysis.  Within  the  anal¬ 
ysis  the  POE  values  are  used  to  omit  errant  data  and  are 
also  used  in  assigning  weights  to  the  data.  The  output  of 
the  analysis  is  given  to  an  ocean  model  either  as  3D  fields 
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Automated  QC  with  Condition  Flags 
(Probability  of  Error) 


Fig.  1.  RNCOM/NCODA  data  flowchart. 


or  3D  increments.  The  ocean  model  assimilates  these 
fields  and  produces  forecast(s)  that  can  then  be  used  for 
future  quality  control  tests  or  for  other  purposes  such  as 
adaptive  sampling.  Additional  details  of  this  flowchart 
will  be  explained  below. 

B.  NCODA  Runtimes 

Part  of  the  NCODA  implementation  process  was  to 
determine  how  many  processors  on  the  Navy  DSRC 
were  optimal  for  running  NCODA  [3].  On  the  DSRC 
system,  each  processor  consisted  of  8  CPUs.  Tests  were 
done  for  three  CPU  configurations  to  determine  how 
the  runtimes  responded.  The  NCODA  code  is  a  mixture 
of  serial  and  parallel  code.  As  such  there  is  a  limit  to 
how  fast  the  code  can  be  run.  In  figure  2  an  example 
is  shown  of  the  runtimes  for  the  NCODA-SA  Global3D 
region  with  8  processors.  Runtimes  are  affected  by  (1)  the 
size  of  state  vector  (proportional  to  the  number  of  grid 
points),  (2)  the  number  of  CPUs,  and  (3)  the  number  of 
input  observations. 


Fig.  2.  NCODA  Analysis  Global3D  8  processors. 

Each  NCODA  region  runs  a  combination  of  2D  and 
3D  serial  and  parallel  processes.  The  preprocessing  and 
postprocessing  is  serial,  while  the  analysis  has  some 


parallel  code.  The  2D  processes  involve  2D  OI  of  SST 
and  SSHA.  The  3D  processes  are  concerned  with  the 
3D  MVOI  analysis  and  its  results.  In  figure  2  the  2D 
preparation  is  shown  in  yellow,  2D  analysis  in  cyan,  2D 
postprocessing  in  magenta,  3D  preparation  in  blue,  3D 
analysis  in  green,  and  3D  postprocessing  in  red.  It  is 
readily  apparent  that  the  majority  of  the  runtime  is  spent 
in  the  3D  MVOI. 

The  following  table  summarizes  the  results  for  8, 
16,  and  32  processors  (64,  128,  256  CPUs).  Shown  is 
the  average  analysis  runtime,  not  including  pre/post¬ 
processing  (the  pre/post-processing  is  serial,  while  the 
analysis  is  partly  parallelized). 


3D  Global  Analysis  Runtime  Improvement 


CPU 

Runtime  (minutes) 

Decrease 

(minutes) 

64 

40 

- 

128 

25 

15 

256 

20 

5 

The  speedup  of  the  NCODA  OI/MVOI  analysis  using 
multiple  processors  in  parallel  is  limited  by  the  time 
needed  for  executing  the  serial  portion  of  the  program 
and  how  efficient  the  parallelizing  is.  Amdahl's  law  [4] 
can  be  used  to  find  the  maximum  expected  improvement 
to  the  NCODA  runtime  when  only  part  of  the  system  can 
be  improved. 

Write  Amdahl's  Law  as  SPEEDUP  —  jzrp,  where 
SPEEDUP  is  the  potential  program  speedup  if  a  P 
fraction  of  code  can  be  parallelized.  The  results  for  the 
runtimes  of  Global2D  OI  and  Global3D  MVOI  are  shown 
in  the  following  table.  Gains  for  the  second  doubling 
are  half  of  the  first,  perhaps  because  (1)  the  parallelizing 
could  be  improved  and  (2)  the  scalability  of  parallelizing 
is  reaching  its  limit.  Based  upon  the  results  shown 
in  the  these  two  tables  and  the  current  configuration 
limitations,  generally  the  NCODA  2D  is  run  with  32  or 
64  CPUs,  while  the  NCODA  3D  is  run  with  128  CPUs. 


REGION 

A  CPU 

Speedup 

Paralleliza 

Global2D 

64  to  128 

1.75 

43  % 

Global2D 

128  to  256 

1.28 

23  % 

Global3D 

64  to  128 

1.62 

38  % 

Global3D 

128  to  256 

1.26 

21  % 

C.  NCODA  OcnQC 

Quality  control  (QC)  of  oceanic  observations  is  an 
important  requirement  for  the  running  of  the  NCODA 
system  [5].  Errant  data  must  be  flagged,  similar  data 
must  be  pooled,  overly  abundant  data  must  be  thinned, 
and  duplicate  data  must  be  eliminated.  In  this  process, 
valid  but  extreme  data  values  are  hopefully  retained.  The 
ultimate  goal  of  the  QC  is  to  feed  "clean"  data  to  the 
analysis  algorithm  so  that  the  output  of  the  analysis  will 
not  be  misleading  and  the  possibility  of  making  wrong 
decisions  is  reduced  or  eliminated. 


QC  of  profile  data  is  of  particular  importance  due  to 
the  need  for  information  at  depth  and  the  small  number 
of  profiles  relative  to  other  data  types.  A  sequence  of 
gross-error  data  checks  are  performed  including  valid- 
value  range  tests,  land-sea  checks  (being  implemented), 
and  location  (speed)  tests  (being  implemented).  A  series 
of  instrumentation  error  checks  are  then  performed, 
such  as  a  sensor  drift  check.  Cross  validation  checks 
are  performed  to  ensure  the  consistency  of  observations 
within  and  between  analysis  variables.  In  the  within- 
variable  consistency  check,  an  OI  analysis  is  performed 
at  the  profile  location,  based  upon  nearby  valid  data  (and 
excluding  the  datum  being  checked). 

The  last  (and  most  important)  check  before  the  anal¬ 
ysis  is  the  background-field  checks.  The  background 
field  checks  can  involve  climatology,  global  and  regional 
analyses,  or  short-term  forecasts.  The  background  fields 
and  background  error  fields  closest  in  time  to  the  datum 
are  interpolated  to  its  position.  The  probability  of  an 
erroneous  value  (the  "probability  of  error"  or  POE)  is 
calculated  from  the  anomalies  of  the  background  fields 
with  respect  to  the  datum  [5]. 

Figure  3  displays  several  of  the  background  fields  used 
in  computing  the  POE,  which  is  computed  both  as  an 
overall  value  and  a  function  of  depth. 

A  final  QC  check  is  done  within  the  analysis  to  re¬ 
move  any  data  that  have  passed  the  previous  checks. 
The  normalized  innovation  corresponding  to  the  data 
is  tested  against  the  analysis  background,  and  if  it  is 
beyond  a  predefined  number  of  standard  deviations, 
the  data  (innovation)  is  rejected.  This  assumes  that  the 
background  error  covariance  (P;,)  and  observation  error 
covariance  (R)  are  reasonably  known.  It  is  best  to  use 
a  high  tolerance  value,  and  four  standard  deviations  is 
commonly  specified. 

D.  NCODA  Analysis 

The  NCODA  3D  MVOI  algorithm  is  formulated  in 
observation  space  as  follows  (with  definitions  in  the 
following  table): 

Xfl  =  Xfa  +  PfcHT  [HP(,Ht  +  R]  -1  [y  -  H(xb)} 


y 

Observation  Vector 

Xfa 

Background  Vector 

x« 

Analysis  Vector 

H,  H 

Forward  Operator,  Matrix1 

R 

Observation  Error  Covariance 

Pb 

Background  Error  Covariance 

y  -H(xb) 

Innovation  Vector 

y-  H(xfl) 

Residual  Vector 

X«  -xfc 

Increment  Vector 

1  The  Forward  Operator  H  is  a  method  of  converting  a 
forecast  model  variable  to  an  observed  variable.  For  NCODA 
the  observations  and  forecast  variables  are  the  same,  and  H 


Fig.  3.  NCODA  OcnQC  profile  plot.  Displayed  are  the  profile  and  the 
anomaly  with  respect  to  the  profile.  The  "probability  of  error"  (POE) 
is  given  both  as  a  function  of  depth  (yellow  curve)  and  as  an  overall 
value. 

reduces  to  a  spatial  interpolation. 

The  analysis  increment  is  equivalently  defined 
by  P/,Ht  [HP,,Ht  +  R]-1[y  —  H  (x/, )].  The  quantity 
PfohfjHPfcHT  +  R]-1  is  the  weight  matrix  (also  commonly 
called  the  Kalman  gain  matrix).  Within  the  3D  MVOI 
algorithm,  the  variables  are  constrained  such  that 
hydrostatic  balance  and  geostrophy  are  maintained. 

The  solution  of  the  MVOI  equation  is  carried  out 
by  an  overlapping  volume  approach  [6].  Eight  volume 
solutions  are  computed  for  each  analysis  grid  point, 
with  each  one  weighted  by  its  distance  from  the  volume 
center.  The  volume  sizes  are  a  function  of  the  local  cor¬ 
relation  length-scale,  such  that  eight  correlation  length- 
scales  are  used.  Smooth  analysis  increments  result  from 
this  combinition  of  overlapping  volumes  and  multiple 
correlation  length-scales  within  a  volume. 

The  background  and  observation  error  covariances 
(Pj,  and  R)  are  separated  into  an  error  variance  and 
a  correlation  [1],  The  correlations  are  further  separated 
into  horizontal  (C/J  and  vertical  (Cs)  components.  A 
tunable  flow-dependent  correlation  is  also  introduced  at 
this  point  (C  f  ).  The  total  background  error  correlation  is 
then  represented  as 

Q,  =  Q,CvCf 

The  default  horizontal  correlation  length-scale  used 
by  NCODA  is  the  (location-dependent)  first  order  baro- 
clinic  Rossby  radius  of  deformation  (see  figure  4).  This 
horizontal  correlation  length-scale  can  be  either  scaled 


or  replaced  by  user-preferred  values  (such  as  those 
produced  by  the  ensemble  method  in  section  IV). 

The  vertical  correlation  length-scales  can  (1)  be  con¬ 
stant,  (2)  monotonically  increase  or  decrease  with  depth, 
or  (3)  vary  with  the  background  vertical  density  gra¬ 
dient.  The  third  option  allows  the  vertical  correlation 
length-scale  to  be  large  when  the  water  column  stratifi¬ 
cation  is  weak,  or  small  when  the  stratification  is  strong. 
It  may  be  represented  as 

hv  —  ps/  {dp/dz) 

where  hv  is  the  resulting  vertical  correlation  length-scale, 
ps  is  a  "change  in  density  stability  criterion"  (which 
defines  a  well-mixed  layer),  and  dp/dz  is  the  vertical 
density  gradient. 

The  flow-dependent  correlation  is  computed  from  a 
scaling  of  the  geopotential  height  difference  between 
two  locations.  This  flow-dependence  affects  both  the 
horizontal  and  vertical  correlations.  Theoretically  other 
fields  that  give  a  good  indication  of  the  flow-field  could 
also  be  used. 

The  NCODA  covariances  can  thus  be  "tuned"  for 
regional  applications  by  using  regionally  dependent  (1) 
horizontal  and  vertical  correlation  length-scales,  and  (2) 
flow-dependence. 


Fig.  4.  NCODA's  location-dependent  default  horizontal  correlation 
length-scale  based  on  the  first  order  baroclinic  Rossby  radius  of 
deformation  (figure  courtesy  of  Jim  Cummings). 

An  important  part  of  the  NCODA  system  is  keeping 
track  of  the  data  and  analysis  mismatches  to  aid  in 
computing  a  reliable  estimate  of  the  background  errors. 
To  facilitate  this,  there  are  three  items  of  particular 
interest  for  which  NCODA  keeps  a  time  history:  the 
data  error,  the  predicted  background  error,  and  an  innovation 
error  check  (the  innovation  error  check  will  be  covered 
in  section  III- A). 

The  data  error  is  a  weighted  running  sum  of  squares 
from  a  time  history  of  the  analyzed  increment  fields. 
These  are  used  to  compute  the  background  errors  of  the 
appropriate  analysis  variables. 

The  predicted  background  error  is  computed  from  a 
weighted  time  history  of  the  analyzed  increment  fields, 
the  data  error  fields,  and  the  climatological  or  model 


error  field.  It  is  also  referred  to  as  a  prediction  or 
forecast  error,  since  it  is  the  error  attached  to  the  MVOI 
background  field,  which  is  commonly  obtained  from  the 
forecast  of  an  ocean  model.  The  predicted  background 
error  (f/, )  is  computed  by 


n 

4  =P  ■  ((E  Wk(Xa  ~  Xb)l)  +  ZVn+ i((xfl  -  Xb )2}) 

k=  1 

+  (!  -  P)  ■  <4 

with  /)  =  exp(— t/tc)2  .  r  is  the  age  of  the  data. 
tc  adjusts  the  rate  at  which  the  background  error  ap¬ 
proaches  the  climatological  or  model  error  in  the  absence 
of  observations.  Time  scales  range  from  wlO  days  for 
surface-  and  mixed-layer  variables  to  «30  days  for 
variables  at  depth,  n  is  the  number  of  days  in  the  past 
to  use  in  computing  the  weighted  summation  (with  k 
being  the  day  index),  are  weight  functions  defined 
by  wb  =  (1  —  0)*-1,  with  cf)  an  adjustable  tuning  factor. 
The  brackets  (...)  represent  a  long-term  mean  increment 
vector.  <7/,  is  the  expected  error  in  the  analysis  variable 
in  the  long-term  absence  of  observations. 

NCODA  keeps  track  of  the  observation  age  (at  all 
model  grid  points)  used  in  computing  the  predicted 
background  error.  An  update  to  any  variable  in  the 
MVOI  analysis  (from  an  observation)  results  in  an  up¬ 
date  to  the  observation  age  of  all  related  variables.  Figure 
5  is  an  example  of  the  observation  age  at  200  meter  in 
the  North  Atlantic.  Since  updates  at  depth  occur  only 
from  in  situ  profiles  or  synthetic  bathythermographs, 
grid  points  at  depth  are  updated  less  often  than  those 
near  the  surface.  However,  there  is  also  less  oceanic 
activity  at  depth,  and  so  the  fields  at  depth  are  closer  to 
climatology  than  those  at  the  surface.  An  example  of  the 
resulting  predicted  background  error  for  SST  is  shown 
in  figure  6,  with  the  largest  errors  occurring  mostly  in 
the  regions  of  largest  mesoscale  activity. 

Prior  to  the  3D  MVOI,  2D  OI  analyses  of  sea-surface 
temperature  (SST),  sea-surface  height  anomaly  (SSHA), 
and  sea  ice  are  performed.  The  2D  OI  preprocessing  of 
the  SST  data  is  done  to  reduce  the  runtime  of  the  3D 
MVOI.  The  2D  OI  of  the  SSHA  data  is  done  in  order 
to  determine  at  which  locations  to  produce  synthetic 
bathythermographs  as  a  means  of  communicating  the 
ocean  topographic  information  to  the  subsurface. 

Before  the  SST  data  are  passed  to  the  3D  MVOI 
(upper-left  plot  in  figure  7),  they  are  preprocessed  by 
a  2D  OI.  The  output  of  the  2D  OI  is  an  SST  increment 
(lower-left  plot  in  figure  7)  to  which  a  threshold  is 
applied  and  then  a  subsampling  is  performed  (upper- 
right  plot  in  figure  7).  The  resulting  SST  data  are  then 
passed  on  to  the  3D  MVOI.  The  resulting  temperature 
increment  at  the  surface  (lower-right  plot  in  figure  7)  is 
seen  to  be  almost  identical  to  the  SST  increment  in  which 


Sampled  SST  (14,228) 


all  the  SST  data  were  used.  This  preprocessing  results  in 
a  very  large  time  savings  in  running  the  3D  MVOI. 

An  important  part  of  NCODA  is  the  transfer  of  oceanic 
topographic  information  to  the  subsurface.  This  is  done 
by  means  of  synthetic  bathythermographic  (BT)  pro¬ 
files.  These  profiles  can  be  created  by  two  methods: 
Cooper/Haines  [7],  and  the  Modular  Ocean  Data  Assim¬ 
ilation  System  (MOD AS)  [8]. 


Observation  Age  (hr)  200  M  Depth 
24  Aug  09  00Z  Tau  000  11  km  grid 


0  40  80  120  160  200  2  10 


Fig.  5.  An  example  of  the  "observation  age"  variable  (at  200-meters) 
which  NCODA  uses  to  determine  how  much  to  relax  a  variable  to 
climatological  values.  The  observation  age  (at  a  grid  point)  is  affected 
by  updates  to  any  of  the  related  MVOI  variables. 
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Fig.  6.  An  example  of  the  "predicted  background  error"  for  SST. 
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Selection  Criteria:  ST  0.1  °C,  sample  every  4th  grid  node 

Fig.  7.  Prior  to  the  3D  MVOI  of  the  NCODA  Analysis,  the  original  SST 
data  are  passed  through  a  2D  OI  that  thresholds  and  subsamples  it 
(figure  courtesy  of  Jim  Cummings). 


The  Cooper/Haines  method  ("direct  method")  needs 
an  ocean  forecast  model.  It  adjusts  the  forecast  den¬ 
sity  profile  to  be  consistent  with  the  change  in  SSHA 
(as  measured  by  the  altimeters).  The  temperature  and 
salinity  are  computed  simultaneously,  and  the  water 
mass  properties  on  subsurface  isopycnals  are  conserved. 
Observation  errors  are  computed  from  prior  estimates 
and/ or  the  analysis  error  variance  plus  the  residual  error 
from  an  iterative  fit  of  the  density  adjustment. 

The  MODAS  method  is  independent  of  an  ocean 
forecast  model.  It  computes  temperature  at  depth 
from  stored  regressions  (empirical  orthogonal  functions 
(EOFs))  between  anomalies  of  climatological  tempera¬ 
ture  and  dynamic  height.  Salinity  values  are  computed 
from  climatological  relationships  between  temperature 
and  salinity.  Observation  errors  are  computed  from  prior 
estimates  and/or  the  regression  residuals. 

Note  that  both  methods  generate  temperature  and  sa¬ 
linity  at  depth  using  SSHA  and  SST  predictor  variables. 
A  strong  point  of  the  Cooper  /Haines  method  is  that 
it  does  not  introduce  spurious  water  masses  into  the 
model.  However,  it  cannot  correct  for  long-term  drift  of 
water  mass  characteristics  in  the  model.  Given  a  good 
ocean  forecast  model,  the  Cooper/Haines  method  may 
be  preferred  [1].  The  RNCOM/NCODA  system  uses  the 
MODAS  method  at  this  time.  When  NCODA  is  run 
in  stand-alone  mode,  the  MODAS  method  is  the  only 
applicable  option. 

Both  methods  make  synthetic  BTs  at  locations  deter¬ 
mined  by  a  2D  OI  of  SSHA  prior  to  the  3D  MVOI.  Figure 
8  shows  an  example  of  an  SSHA  increment.  At  locations 
where  the  absolute  value  of  the  increment  is  larger  than 
a  preset  "noise"  value  (for  example  2.0  centimeters), 
a  cross-hatched  subsampling  will  be  performed  (see 
figure  9).  The  resulting  grid  spacing  is  also  adjustable. 
By  means  of  these  two  parameters  (along  with  the 
correlation  length-scales)  the  influence  of  the  altimetry 
can  be  adjusted. 


II.  NCODA  Coupled  with  GNCOM 
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The  Global  Navy  Coastal  Ocean  Model  (GNCOM)  [9] 
makes  use  of  several  ocean  models  and  data  assimilation 
systems  to  accomplish  its  prediction  of  ocean  forecasts. 
These  include  the  Global  Navy  Layered  Ocean  Model 
(GNLOM)  [10],  the  2D  OI  capability  of  MODAS  [11], 
the  3D  OI  capability  of  MODAS  [8],  and  the  MVOI 
method  of  NCODA  [1].  The  relationship  of  these  systems 
is  shown  in  figure  10. 

GNCOM  encompasses  the  open  ocean  to  5  m  depth 
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Fig.  8.  An  example  of  a  2D  SSHA  increment  output  from  the  2D  OI 
preprocessing  of  the  NCODA  Analysis  . 
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Fig.  9.  After  applying  a  height  threshold  test  to  the  2D  SSHA  increment, 
a  grid  of  locations  is  generated,  which  will  dictate  where  synthetic  BTs 
will  be  created  for  input  to  the  3D  MVOI  analysis. 


Fig.  10.  GNCOM/NCODA  data  flowchart. 


in  a  curvilinear  global  grid  with  1/8  degree  grid  spacing 
at  45  N,  extending  from  80  S  to  a  complete  arctic  cap 
with  grid  singularities  mapped  into  Canada  and  Russia 

[12] ,  The  model  employs  40  vertical  Sigma/Z  levels,  with 
Sigma  in  the  upper  ocean  and  coastal  regions,  and  Z  in 
the  deeper  ocean.  The  real-time  systems  (GNCOM  and 
GNLOM)  are  forced  with  the  global  1/2°  Navy  Opera¬ 
tional  Global  Atmospheric  Prediction  System  (NOGAPS) 

[13] . 

The  input  observations  to  GNCOM  follow  two  dif¬ 
ferent  processing  paths,  with  the  result  of  each  path 
producing  3D  temperature  and  salinity  profiles.  Satellite 
SST  are  processed  by  the  2D  OI  of  MODAS  2D.  Altimetry 
is  input  to  GNLOM,  which  outputs  2D  fields  of  SSHA 
and  mixed  layer  depth  (MLD).  These  are  combined  with 
the  2D  SST  output  of  MODAS  2D  to  form  the  inputs 
to  MODAS  3D.  MODAS  3D  then  creates  global  3D 
synthetic  temperature  and  salinity  fields  that  become  the 
background  (first  guess)  field  of  the  NCODA  3D  MVOI. 

The  in  situ  temperature  and  salinity  observations  from 
gliders,  buoys,  etc.  are  processed  by  the  NCODA  OcnQC 
(section  I-C).  These  correspond  to  the  observation  vector 
mentioned  previously  (y  in  the  3D  MVOI  equation). 
The  background  field  produced  by  MODAS  3D  (H(x j,)) 
is  subtracted  from  these  observations  to  produce  the 
innovations  that  are  transformed  by  the  3D  MVOI  into 
the  output  increment.  These  innovations  are  computed 
without  the  First  Guess  at  Appropriate  Time  (FGAT) 
method  of  section  III.  This  increment  is  added  to  the 
input  background  field  to  produce  a  3D  temperature  and 
salinity  field.  GNCOM  is  relaxed  to  this  field  over  a  three 
day  hindcast. 

This  global  ocean  prediction  strategy  [14]  satisfies  the 
need  for  ocean  models  to  have  a  high  horizontal  reso¬ 
lution  (required  to  successfully  simulate  the  variability 
of  mesoscale  features)  and  a  high  vertical  resolution 
(required  near  the  surface  to  resolve  the  physics  of  the 
upper  ocean)  [15].  This  first  generation  global  ocean  pre- 


diction  system  meets  these  needs  by  combining  GNLOM 
(for  high  horizontal  resolution)  and  GNCOM  (for  high 
vertical  resolution). 

GNCOM  produces  forecasts  of  up  to  72  hours.  Its 
output  is  used  as  boundary  and  initial  conditions  for 
nesting  the  higher  resolution  RNCOM/NCODA  system 
(to  be  discussed  in  section  III)  inside  of  it.  Since  GNCOM 
does  not  include  tides,  these  are  introduced  from  an 
external  tidal  database  when  GNCOM  is  coupled  with 
a  nested  RNCOM/NCODA  system. 
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Fig.  11.  The  RNCOM/NCODA  system  uses  an  "incremental  (analysis) 
updating"  (or  data  insertion)  method  of  assimilating  the  temperature 
(T)  and  salinity  (S)  increments.  The  T,S  increments  are  gradually 
inserted  during  a  24-hour  hindcast.  This  helps  maintain  dynamic 
balance  and  cut  down  on  spurious  gravity  waves. 


III.  NCODA  Coupled  with  RNCOM 

The  Regional  Navy  Coastal  Ocean  Model  (RNCOM) 
[16]  is  the  Navy's  choice  for  high  resolution  applications. 
It  is  based  on  the  Princeton  Ocean  Model  [17],  combining 
Sigma-layers  and  Z-levels  [18],  permitting  use  of  a 
hybrid  vertical  coordinate  system  [15].  The  vertical  grid 
is  set  up  to  offer  a  choice  of  Sigma-layers  or  Z-levels, 
or  some  combination  with  Sigma-layers  in  the  shallow 
water  and  Z-levels  in  the  deeper  water.  A  typical  imple¬ 
mentation  of  RNCOM  may  be  at  1/32  or  1/64  degree 
resolution,  with  45  levels  (15  Sigma-layers  and  30  Z- 
levels),  and  forecasts  from  72  to  96  hours.  The  model  is 
forced  with  the  high  resolution  regional  Coupled  Ocean 
Atmosphere  Mesoscale  Prediction  System  (COAMPS) 
[19]. 

RNCOM  takes  its  boundary  conditions  (BC)  from 
GNCOM.  Since  GNCOM  does  not  include  tides,  they  are 
introduced  by  adding  tidal  heights  and  transports  to  the 
BC  (from  the  Oregon  State  University  tide  model  [20]). 
Internal  to  an  RNCOM  domain,  tidal  forcing  is  done  via 
the  tidal  potential. 

NCODA  is  the  data  assimilation  system  used  by  RN¬ 
COM.  The  data  flow  of  the  RNCOM/NCODA  system 
was  shown  in  figure  1  of  section  I-A.  All  of  the  raw 
observations  pass  through  the  NCODA  OcnQC  system 
as  detailed  in  section  I-C.  With  RNCOM  the  innovations 
are  computed  using  the  forecast  from  yesterday  that 
is  closest  in  time  to  the  observation  (the  First  Guess 
at  Appropriate  Time  (FGAT)  method).  The  background 
(first  guess)  field  of  the  MVOI  is  the  24-hour  forecast 
from  yesterday.  The  temperature  and  salinity  increment 
fields  produced  by  the  3D  MVOI  are  what  is  used  by 
RNCOM  (in  contrast  to  GNCOM). 

In  order  to  minimize  spurious  gravity  waves  (result¬ 
ing  from  introducing  an  unbalanced  increment  into  the 
model),  the  RNCOM/NCODA  system  uses  an  "incre¬ 
mental  (analysis)  updating"  (or  data  insertion)  method 
of  assimilating  the  temperature  and  salinity  increments 
(see  figure  11).  This  is  accomplished  by  a  24-hour  hind- 
cast  that  starts  with  yesterday's  nowcast  and  gradually 
inserts  the  temperature  and  salinity  increments.  The 
velocity  fields  are  allowed  to  adjust  via  geostrophy. 


A.  RNCOM/NCODA  Metrics 

The  performance  of  the  RNCOM/NCODA  system  can 
be  monitored  by  several  graphics.  In  figure  12  are  shown 
profile  plots  of  the  temperature  fields  from  the  in  situ 
observation  (Obs),  NCODA  analysis  (Anal),  RNCOM 
nowcast  (NCST)  for  today,  RNCOM  24-hour  forecast 
(FCST)  from  yesterday,  GDEM  climatology  (Clim),  and 
the  background  field  (BG)  with  its  associated  error  (E). 
Recall  that  the  NCODA  analysis  increment  is  passed 
to  the  RNCOM  model,  which  gradually  inserts  it  over 
a  24-hour  hindcast  to  produce  the  nowcast  for  today. 
The  background  field  (BG)  is  the  FGAT  field  used  when 
computing  the  innovation  for  input  to  the  analysis. 
Figure  13  is  an  anomaly  plot  of  the  same  quantities 
with  respect  to  the  in  situ  observation.  Note  that  in 
the  anomaly  plot  the  "background  anomaly"  is  actually 
the  "innovation",  which  is  used  by  the  analysis  with 
appropriate  weighting  based  on  its  probability  of  error 
values. 
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Fig.  12.  Profile  plots  of  the  temperature  fields  from  the  in  situ  obser¬ 
vation,  NCODA  analysis,  and  the  RNCOM  nowcast  and  forecast. 
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Fig.  15.  The  residual  RMS  of  the  NCODA  temperature  analysis  over 
time  (as  a  regional  average)  is  studied.  The  residual  RMS  is  largest 
between  100  and  300  meters  in  this  region  around  Hawaii,  with 
maximum  variations  likely  occurring  where  internal  tides  have  the 
most  effect. 


Fig.  13.  Profile  plots  of  the  anomaly  fields  (difference  from  the  in 
situ  observation)  from  (1)  the  NCODA  analysis  and  (2)  the  RNCOM 
nowcast  and  forecast.  Overall  RMS  values  are  shown  within  the 
legend.  Note  that  the  background  anomaly  (B-O)  is  the  innovation 
computed  using  FGAT.  In  general  the  analysis  should  have  the  best  fit 
to  the  data,  but  the  nowcast  should  have  the  best  dynamical  balance. 

How  well  the  system  is  removing  any  bias  and  re¬ 
ducing  the  root-mean-square  (RMS)  variability  is  also 
monitored.  Scatter  plots  are  made  showing  the  results 
before  and  after  the  analysis,  which  correspond  to  the 
statistics  of  the  innovations  (before)  and  the  statistics  of 
the  residuals  (after).  Figure  14  is  an  example  of  a  regional 
SST  analysis  where  a  large  bias  was  corrected  and  the 
RMS  was  significantly  reduced.  An  example  of  tracking 
the  regional  RMS-at-depth  over  time  is  shown  in  figure 
15.  In  this  case  for  the  Hawaii  region,  the  RMS  signals 
seem  to  contain  a  48-hour  signal,  with  the  largest  RMS 
values  occurring  at  depths  where  the  internal  tides  have 
the  most  effect. 
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Fig.  14.  Scatter  plots  of  matchup  statistics  before  (blue)  and  after  (red) 
the  NCODA  Analysis. 

IV.  RNCOM /NCODA  Ensembles  and  ASAP 

Recall  in  section  I-D  that  NCODA  is  able  to  introduce 
flow-dependence  into  the  data  assimilation  by  means  of 
a  flow-dependent  correlation  (C  f  ).  An  improved  method 


of  determining  the  flow-dependence  is  to  have  an  en¬ 
semble  of  ocean  model  runs  valid  at  the  same  times. 
This  will  allow  a  better  estimation  of  the  background 
error  (variance). 

The  results  of  ensemble  forecast  runs  will  not  only 
give  better  estimates  of  background  errors,  but  can  also 
provide  estimates  of  future  errors  (in  sound  speed  vari¬ 
ability,  etc.),  which  can  be  used  to  predict  where  obser¬ 
vations  should  be  taken  now  in  order  to  reduce  those 
errors.  This  is  important  when  the  Navy  must  decide 
how  and  when  to  deploy  a  limited  number  of  assets 
in  order  to  attain  a  better  knowledge  of  the  battlespace 
environment. 

To  accomplish  this,  the  RNCOM/NCODA  system  is 
being  used  as  part  of  an  Adaptive  Sampling  and  Pre¬ 
diction  (ASAP)  project  that  makes  use  of  the  Ensemble 
Transform  (ET)  and  Ensemble  Transform  Kalman  Filter 
(ET  KF)  applied  to  ensemble  runs  of  the  RNCOM. 

Letting  e  represent  the  error  (uncertainty)  in  a  quan¬ 
tity,  the  error  in  a  RNCOM  forecast  may  be  represented 
as  follows: 

Not  Yet  Included 

, - ~ - c 

£ Forecast  =  £ Forcing  "F  OC  ~F  €BC  ~F  £  Model  ~F  e Turbulence 

There  are  other  terms  that  could  be  included,  such  as  er¬ 
rors  in  the  bathymetry  and  sub-grid  variability.  Only  the 
first  two  sources  of  error  are  implemented  at  this  time. 
Thus,  the  ensemble  generation  for  RNCOM/NCODA  is 
done  by  perturbing  the  forcing  and  initial  conditions. 

Perturbation  of  the  forcing  takes  the  form  of  a  space- 
time  deformation  of  the  atmospheric  forcing  [21].  This 
results  in  an  ensemble  of  atmospheric  states  -  an  inde¬ 
pendent  atmospheric  forcing  for  each  RNCOM  ensemble 
member. 

Perturbation  of  the  initial  conditions  (IC)  is  done  by 
the  Ensemble  Transform  (ET)  technique  [22]  [23].  The 
ET  method  uses  the  best  available  estimate  of  analysis 


error  covariance  to  transform  forecast  perturbations  into 
analysis  perturbations  by  finding  K  distinct  linear  com¬ 
binations  of  K  forecast  perturbations  that  (1)  are  equally 
likely,  (2)  lie  within  the  vector  subspace  of  forecast  per¬ 
turbations,  (3)  are  quasi-orthogonal  (although  they  sum 
to  zero),  and  (4)  have  expected  squared  amplitudes  equal 
to  the  trace  of  the  best  available  estimate  of  the  analysis 
error  covariance  matrix.  For  the  current  work,  the  best 
available  estimate  of  the  analysis  error  covariance  is 
obtained  from  the  NCODA  analysis  error. 

The  K  perturbations  to  the  IC  are  then  added  to  the 
RNCOM  control  run  to  produce  K  initial  states.  These  K 
initial  states  (ensemble  members)  are  then  integrated  for¬ 
ward  in  time  with  their  respective  atmospheric  forcings 
to  produce  K  forecasts.  These  forecasts  are  used  with  the 
ET  KF  technique  for  ASAP  purposes  [24].  At  this  point 
the  process  can  be  restarted  for  the  next  ASAP  run. 

The  use  of  the  RNCOM/NCODA  system  for  ASAP 
has  been  done  for  several  ocean  measurement  exercises 
[25]  [26]  [27]  [28].  To  make  the  ability  to  perform  ASAP 
and  adjust  the  parameters  it  uses  more  accessible  to 
the  warfighter,  a  software  system  called  "Targeted  Ob¬ 
servations  for  Forecast  Uncertainty  (TOFU)"  has  been 
developed  [25]  [27].  It  creates  metrics  that  can  be  used 
for  ASAP,  and  that  also  assess  how  well  the  ASAP  is 
performing. 

The  TOFU  system  produces  "Target  Observations 
Summary  (TOSum)"  maps  as  shown  in  figure  16.  These 
display  the  relative  impact  of  each  possible  temperature- 
salinity  profile  observation  in  reducing  the  predicted  er¬ 
rors  of  a  set  of  target  variables  (such  as  sonic  layer  depth 
and  below  layer  gradient)  over  a  user-specified  target 
area  and  for  a  given  target  forecast  time.  These  maps 
can  then  be  integrated  with  other  criteria  to  optimally 
design,  deploy,  and  update  local  observation  networks 
that  will  provide  the  best  accuracy  of  ocean  dynamic 
inputs  into  end  user  applications. 

In  the  case  of  figure  16,  the  red  box  outlines  the  target 
area.  There  are  four  gliders,  with  three  within  the  target 
area.  The  possible  deployments  of  each  glider  (taking 
account  of  ocean  currents)  are  shown  by  white  lines.  Red 
areas  on  the  map  indicate  locations  where  measurements 
will  have  a  large  relative  impact  on  reducing  the  forecast 
error  (uncertainty)  within  the  target  region.  Blue  areas 
will  have  a  small  impact. 

V.  Future  NCODA  Enhancements 

The  data  assimilation  algorithm  used  by  NCODA  is 
being  updated  in  several  stages.  The  first  transition 
will  be  to  a  3D  variational  (3DVar)  approach  [29].  The 
equation  to  be  solved  is  the  same  as  the  one  given 
in  section  I-D,  but  using  variational  techniques.  The 
techniques  will  come  from  the  NRL  Atmospheric  Vari¬ 
ational  Data  Assimilation  System  (NAVDAS),  which  is 
an  observation-space  based  3DVar  system  for  generating 
atmospheric  state  estimates.  NAVDAS  permits  the  direct 
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Fig.  16.  TOSum  map  for  three  gliders,  showing  possible  deployments 
of  the  gliders  and  the  relative  impact  of  each  upon  the  target  area  (red 
box). 

assimilation  of  satellite  radiances,  which  have  resulted 
in  significant  forecast  improvements  at  other  forecast 
centers  [30]. 

The  next  transition  will  be  to  a  4D  variational  (4DVar) 
approach  [31]  applied  to  the  RNCOM  model  [32].  This 
4DVar  approach  will  be  based  upon  the  "Cycling  Repre¬ 
senter  Method".  In  this  method  RNCOM  is  linearized 
about  a  background  state  using  tangent  linearization. 
The  stability  of  this  tangent  linearized  model  (TLM)  is  a 
very  sensitive  function  of  the  background  state,  the  level 
of  nonlinearity  of  the  model,  open  boundary  conditions, 
and  the  complexity  of  the  bathymetry  and  flow  field. 
Based  on  the  TLM  stability  time  period,  the  Representer 
Method  is  cycled  by  splitting  the  time  period  of  the 
assimilation  problem  into  short  intervals.  The  interval 
time  period  needs  to  be  such  that  it  is  short  enough 
for  the  TLM  to  be  stable,  but  long  enough  to  minimize 
the  loss  of  information  due  to  reducing  the  temporal 
correlation  of  the  dynamics  and  data. 
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